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We use linear analysis to examine the effect of cosmic rays (CRs) on the 
Parker- Jeans instability of magnetized self-gravitating gaseous disks. We adopt 
a slab equilibrium model in which the gravity (including self-gravity) is perpen- 
dicular to the mid-plane, the magnetic field lies along the slab. CR is described as 
a fluid and only along magnetic field lines diffusion is considered. The linearised 
equations are solved numerically. The system is susceptible to Parker- Jeans in- 
stability. In general the system is less unstable when the CR diffusion coefficient 
is smaller (i.e., the coupling between the CRs and plasma is stronger). The sys- 
tem is also less unstable if CR pressure is larger. This is a reminiscence of the 
fact that Jeans instability and Parker instability are less unstable when the gas 
pressure is larger (or temperature is higher). Moreover, for large CR diffusion 
coefficient (or small CR pressure), perturbations parallel to the magnetic field 
are more unstable than those perpendicular to it. The other governing factor 
on the growth rate of the perturbations in different directions is the thickness of 
the disk or the strength of the external pressure on the disk. In fact, this is the 
determining factor in some parameter regimes. 



Subject headings: cosmic rays — instabilities — ISM: magnetic fields — MHD 
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1. INTRODUCTION 

In our Galaxy (and probably other galaxies as well), the energy densities of thermal gas, 
magnetic field, cosmic rays and turbulence are of the same order of magnitude (~1 eV cm -3 , 
see e.g., Gaisser 1990). Rough equipartition exists between these components indicates that 
they are tightly coupled together, and all of them should have significant influence on the 
structure and evolution of the interstellar medium (ISM). However, the dynamical role of 
CRs on ISM has not been taken seriously in most ISM related models until rather recently 
(see e.g., Ferriere 2001). 

The ISM is highly inhomogeneous and have many different phases of gas clouds, e.g., 
hot phase, warm phase, diffuse clouds, molecular clouds, etc. This is due partly to the ex- 
istence of many instabilities in ISM conditions. Gravitational instability (Jeans instability) 
and magnetic buoyancy instability (Parker instability, Parker 1966) have been considered 
to be responsible for the formation of some large scale clouds (e.g., HI clouds, molecular 
clouds). As ISM is aggregated, filament-like structures are formed by Parker- Jeans insta- 
bility. Extensive analytical works on Parker- Jeans instability have been carried out in the 
past two to three decades (e.g., Elmegreen & Elmegreen 1978; Elmegreen 1982a,b; Hanawa 
et al. 1992; Nagai et al. 1998; Lee & Hong 2005). Recently, this instability has been studied 
by the magnetohydrodynamic (MHD) simulations (e.g., Chou et al. 2000; Lee et al. 2001; 
Kim et al. 2002). Nagai et al. (1998) showed that the angle of the filament formed by this 
instability with respect to the magnetic field lines depends on the strength of the external gas 
pressure of the initial equilibrium gas layer, and Chou et al. (2000) confirmed their results 
by MHD simulations. Nevertheless, the dynamical effect of CRs has not been addressed in 
these works. 

In his original work on the Parker instability, Parker realised the importance of CRs 
on the instability and did consider their effect in a simplified way (Parker 1966). However, 
studies on the instability with more rigorous treatment of the dynamics of CRs appear only 
rather recently. Both linear analysis (Kim et al. 1997; Kim & Hong 1998; Hanasz & Lesch 
1997; Hanasz 1997; Ryu et al. 2003) and MHD simulation (Hanasz & Lesch 2000, 2003; 
Hanasz et al. 2004a,b; Kuwabara et al. 2004) have been performed on models in which the 
diffusion of CRs through plasma is included. These models are based on the hydrodynamic 
description of CRs (e.g., Drury & Volk 1981; Drury 1983; Ko 1992), which is a fairly good 
approximation in studying the structure and evolution of the system, although the spectrum 
of CR is compromised. As far as we know, the effect of the CR diffusion on the Parker- Jeans 
instability has not been studied yet. 

In this paper, the effect of CRs on the Parker- Jeans instability is examined by linear 
analysis. Hydrodynamic approach to CR is adopted and only along field line diffusion is 



- 3- 



considered. As a model for Galactic gaseous disk or compressed molecular cloud, a self- 
gravitating gas layer or slab threaded by magnetic field perpendicular to the gravity is used 
for the unperturbed equilibrium state. The dispersion relations are obtained numerically. 
The paper is organised as follows. In §2, we present the CR-plasma system and the equilib- 
rium model. § 3 is devoted to linear stability analysis, and the results are described in § 4. 
§ 5 provides a summary and discussion. 



2. MODELS 

2.1. Two-Fluid Model 

We adopt a two- fluid model which comprises cosmic rays and thermal gas (e.g., Drury 
& Volk 1981; Drury 1983; Ko 1992). CR is treated as a massless gas with significant energy 
density or pressure. CRs are coupled to the thermal gas via hydromagnetic irregularities or 
waves embedded in the flow. To a first approximation, the only effect of the irregularities 
or waves is contained in the diffusion coefficient of CR. For simplicity, we consider diffusion 
along magnetic field lines only and ignore the cross-field-line diffusion (see e.g., Giacalone & 
Jokipii 1999). Thus k = Kj bb, where b is the unit vector along the magnetic field. 

Assuming infinite conductivity for the plasma (i.e., ideal MHD), the system is governed 
by the total mass equation, momentum equation, thermal and CR energy equations, 

^ + V .(pV) = 0, (1) 



|(PV) + V 



B 2 \ BB 

pVV + P g + P c + — I 



2/i / /io 

+p [v(v> + ip e ) + 2n x v + n x (n x r)] = o , (2) 

dP B 



dt 

dR 



+ V.VP g + 7 g P g V-V = 0, (3) 



c 



M + V • VP C + 7c PV • V - V • ( K ,|bb • VP C ) = , (4) 

supplemented by the Poisson equation for self-gravity and induction equation for magnetic 
field 

VV = 47rGp , (5) 
— = Vx(VxB), (6) 
V • B = . (7) 
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In these equations p, V, P g and 7 g denote the thermal gas density, flow velocity, pressure and 
polytropic index, respectively; P c , 7 C and k\\ denote the CR pressure, polytropic index and 
diffusion coefficient along magnetic field lines, respectively; ip and ip e denote the potential 
of the self-gravity of the gas and the potential of some external gravitational source, respec- 
tively; B is the magnetic field strength and O is the angular velocity of the disk. In the rest 
of the paper we assume that the centrifugal force is balanced by the external gravitational 
source. 



2.2. Equilibrium Model 

We are interested in the stability of a uniformly rotating gaseous disk or a slab. We 
adopt a local Cartesian coordinate system (x, y, z) such that e x = e^, e y = — e r and e z = 
e z , where (r, 0, z) is the cylindrical coordinate system of the disk. All quantities of the 
unperturbed equilibrium state depend on z only. Moreover, assume initial magnetic field is 
in the azimuthal direction, i.e., B = B(z)e x . For simplicity ignore the external gravitational 
force in the z-direction, the equilibrium model satisfies the magnetohydrostatic equation 

Substitute into the Poisson equation (Eq. (5)), we obtain 



_d 

cU 



%P t dz\ ^C s 2 (l + a + /3) ' {) 



where a and (5 are the ratio of the magnetic pressure and CR pressure to the gas pressure, 
P t = (1 + a + P)Pg is the total pressure, and C s = ("f g P g / p) 1 / 2 is the gas sound speed. 
Moreover, we assume a temperature distribution 

T ( z) = l T *> tor\z\<z d , 
I oo , for \z\ > Zd , 

(Hanawa et al. 1992), where Td and z d are the temperature and the half-thickness of the 
(cold) gaseous disk. Large (small) Zd can be interpreted as small (large) external pres- 
sure. Integrate Eq. (9) with this temperature distribution gives the equilibrium pressure and 
density distributions. The density scale height at z = (mid-plane of the gaseous disk) is 
H = C s0 a/ (1 + oi + /5)/(27rGpo7g) , where the subscript denotes the value at the mid-plane. 
Naturally, we normalise the density, velocity and length of the gaseous disk to po, C s o and 
Hq = C s0 / \j2ixGpo r ) g . The CR diffusion coefficient is normalised to C s0 H . In this paper 
we choose 7 g = 1 and 7 C = 4/3. We pick a — 1, f3 — 1 and z d = H as the fiducial model. 
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Fig. 1 shows the distribution of the equilibrium gas pressure and the gravity in the case of 
a = 1 and (5 — 1. 

Parker- Jeans instability is one of the possible mechanisms for the formation of large 
ISM clouds (e.g., molecular clouds, HI clouds), thus the typical values for p , C s0 , H and 
C s0 H are 1.67 x 1(T 21 kg m~ 3 , 5 km s _1 , 200 pc and 3 x 10 22 m 2 s" 1 . The effective sound 
speed C s oVl + a + P ~ 8-5 km s _1 for a — f3 — 1 (Hanawa et al. 1992). The value of CR 
diffusion coefficient is estimated to be 3 x 10 24 m 2 s -1 (e.g., Brezzinskii et al. 1990; Ptuskin 
2001; Ryu et al. 2003). Thus the fiducial value of normalised k» is 100. On the other hand, 
Parker- Jeans instability may cause the fragmentation and collapse of large ISM clouds. The 
typical values for p , C s o, H and C s qH become 1.67 x 10~ 19 kg m -3 , 5 km s^ 1 , 20 pc and 
3 x 10 21 m 2 s" 1 . In this case the CR diffusion coefficient is smaller, say 3 x 10 23 m 2 s~ x , then 
the fiducial value of normalised k\\ is still 100. 



LINEAR STABILITY ANALYSIS 



3.1. Equations for Perturbations 



We perform standard linear stability analysis on the system of equations Eqs. (l)-(7) 
of the slab equilibrium model described in previous section. Since the unperturbed state 
depends on z only, we assume the perturbed quantities to be of the following form 







Sp 


sv x 




iSV x 


SVy 




i5V y 






sv z 


SP g 




SP g 


SPc 




SPc 


SB X 




5B X 


5By 




5B y 


SB Z 




-i5B z 






5tp 



exp(<rt + ik x x + ik y y) , 



(11) 



where the barred quantities are functions of z only, a is the growth rate, and k x , k y are 
the wave numbers in the x- and ^/-directions, respectively. The set of linearised equations 
becomes 

- kJV X - kySVy + ^ + ~%V Z = , (12) 

p dz pdz 

k - k 1 dB 

<t5V x + —SP, + —5P C + -^5B Z - 2Q5V V + kjip = , (13) 

p p pop dz 
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aSV y + %P g + ^5P C - —5B V + ^-SB X + 2QSV X + kM = , 
P P Pop Pop 



(14) 



- ld5P R ld5P c 1 d ( BbB, 



p p dz p dz \ /j,q 



±±(p + P + * 

p 2 dz 1 s c 2/i 



o/ + dz 



PoP 



0, 



dP - 
a5P g + —^5V Z - 
dz 



l aS p+^SV z 



P 



d;2 



0, 



dP 

a5P c + — 
dz 



7c^c f . dp - \ 2 = reiifca-dP, 



P 



— — —^-5B Z 
B dz 



d 2 # 



+ fcj)<fy + — f - 4ttG5p = , 

^B* + fc XJ B5T4 + ^ W 2 - ^ (a6p + ^<5V 2 ) = , 

obBy + fc x S<jy w = , 
a5B z + k x B5V z = , 
d5B z 

k x 5B x + k y 5B y - = . 



(15) 
(16) 

0, (17) 

(18) 

(19) 

(20) 
(21) 

(22) 



After some manipulations, we obtain the following set of ordinary differential equations, 



d 

dz 







A n 


An 


A 13 







2/1 


2/2 




An 


A22 


A 23 


A24 




2/2 


2/3 













-1 




2/3 


2/4 




An 


A 42 


A 43 







2/4 



(23) 



where the running variables are defined as 



2/1 = 



pSV z 



(24) 



y 2 = 5P t = 5P g + 5P C + —5B X , 

Po 



V3 = Sip, 



2/4 = 5g g 



d5jj 
~~dz~ 



(25) 
(26) 
(27) 



and the matrix elements are 



An = 



e 2 (i + r 2 ) 



x 



1 

D 



a 



+ 



£ 2 £ 2 (i + r 2 ) 



Tr odlnB ldP c 
Vl— — + - 



dz 



r 2 + ^p 

k x 

t2 + t Lt 

fox 



d hip 

dz 



p dz {(1 + K) 



+ V, 



V k 2 x vlv 2 
£ 2 + £2(1 + r 2 ) 



12 £ 2 (i + r 2 ) + J D 



+ 



v^ 2 



x 



£2 

£2 + £2(1 + r 2 ) 



£2(1 + r 2 ) 

K 2 A;2 



(fig + fij)p 

13 ~£ 2 (i + r 2 ) J D£2(i + r 2 ) 




<7 



+ 



1 +? r 
■p2 _ p 



£2 £2(i + r 2 ) 



121 



D 



+ vX 

Tr9 d\nB 

A dz 



ldP r 



C 2 C 



p dz {(1 + K) 



a 

£2 



fcX 2 r 2 
£2(i + r 2 ) 



A 41 = 



AttG 



in which 



V, 



A - 9z 
D 

^23 = 

,dlnS ldP r 



£2 £2(1 + r 2 ) 

P£ 2 (i + r 2 ) 
Ma = p , 



p2 _ K_y 

k x 
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fox. 



dz 
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+ - 



p dz 
~D~ 



(1 + K) 



° 2 | 

£2 + £2(1 + r 2 ) 



£2 + £2(i + r 2 ) 



p2 _ p 

k v 



A 43 = -(k 2 x + k 2 y ) + 



4^ V A 2 A;2p 

d £2(i + r 2 ) 



K 



D = C 2 + 



a 2 + Vlk. 



1 + r r 

kx. 



2aQ 



£2 



Cl 



(1 + K) 



+ 



V 2 k 2 x T 2 



£ 2 £2(1 + r 2 ) 
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£ 



i c p c 



vj 



B 2 



ldP t 

9z = — 



(39) 



P P PoP P dz 

For completeness, we list the perturbed quantities in terms of the running variables (24)-(27), 



sp = -\v£ 



+ 



,dhxB ldP c 



dz 



+ - 



CI 



p dz {(1 + K) 



+ V 2 A 



a 2 
^ + 



a 



+ 



V 2 k 2 



SV, 



5Vy = 



s 2 s 2 (i + r 2 ) 

ak x 

~s 2 (i + r 2 ) 



r 2 _ 



A , 



s 2 1 s 2 (i + r 2 ) 

Vlk 2 x p 



V3 + 



d s 2 (i + r 2 ) 

y-2 



d hip 

dz 



k. 



P 



E 2 (i + r 2 ) 



i + 



V3 + 



m 
p 



, krr 

k,. 



+ V 2 T- X 



1 dp 5p 

— -r-Vi 

fr dz p 



1 dp 5p 

p z dz p 



a 



SV Z = 7/i , 

P 



SP g = C 2 5p - 



5Pc 



a 



c 2 c i 



2 dp 

! dz 

c 2 



dp, 

dz 
dp 



(1 + K) dz 
5g z = 2/4 , 



yi, 

dPc 

dz 



yi , 



SB X = - 



B fdlnB dlnp 



dz 



dz 



yi + B B5V X , 

P ° 



A 



5 B y = ~B8V y , 



5B Z = --BJV Z , 
a 



m_ 

D 



m 

D 



(40) 

(41) 

(42) 
(43) 

(44) 

(45) 
(46) 
(47) 

(48) 
(49) 



By setting up suitable boundary conditions at z = and z = Zd, the set of equations 
(Eq. (23)) is solved for given (k X) k y ). The growth rate is the corresponding eigenvalue that 
satisfies the boundary conditions. 



3.2. Boundary Conditions 

The slab geometry and layered equilibrium model (Eq. (10)) we use in §2 is a kind of 
disk-halo model. We adopt the boundary conditions used by Nagai et al. (1998), Hanawa et 
al. (1992) and Chou et al. (2000) for the interface between disk and halo (at z = z d ). The 
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boundary conditions are derived from the continuity of velocity, pressure, magnetic field, 
gravitational potential, and gravitational acceleration across the interface z = zj, 



V2(Zd) 



9z{z d ) 



kl 



vl 



yi(zd) , 



and 



Vi{z d ) = -^Gyi{z d ) + sjk 2 x + kl y 3 (z d ) 



(50) 



(51) 



The natural boundary conditions at the mid-plane (z = 0) are, 



?/i (0) = 2/4(0) = , (in the case of symmetry) , 
2/2(0) = 2/3(0) = , (in the case of antisymmetry) 



(52) 



4. RESULTS 

In this section we present the dispersion relation of the Parker- Jeans instability with 
respect to different parameters of the system, such as the CR diffusion coefficient k\\, the 
initial ratio of the CR pressure to gas pressure (3, the initial half-thickness of the gaseous 
disk Zd, and the angular velocity Q. We choose a — 1, (5 — 1, Zd/H = 1, k\\/C s qH = 100, 
Q = as the fiducial model (see §2 for the normalisations). The values of z& and k\\ in the 
following discussion are the normalised values. 

We only present the results of symmetric boundary at z — 0, because the growth rate 
in the case of antisymmetric boundary at z = is much smaller (e.g., Chou et al. 2000). 

In the following, a x and o-^max denote the growth rate and maximum growth rate when 
k y = 0, i.e., no perturbations perpendicular to the unperturbed magnetic field; and a y and 
o"j/,max denote the growth rate and maximum growth rate when k x = 0, i.e., no perturbations 
parallel to the unperturbed magnetic field. k Xfir - lt {k yfiT \t) denotes the critical wavenumber 
beyond which the system becomes stable (i.e., a x < 0) when k y = (k x = 0). We are 
particular interested in these two perturbations (k y = and k x = 0), because from the 
numerically studies we have so far, the most unstable mode of the system always occurs at 
either k y = or k x = 0. That means the magnetised gaseous disk tends to break up into 
elongated fragments with axes either perpendicular to {k y = 0) or parallel to (k x = 0) the 
initial magnetic field. 
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4.1. Dependence on CR Diffusion Coefficient km 

Fig. 2 (top) shows the dispersion relation of the case k y = for different k\\ (the CR 
diffusion coefficient). The other parameters take the fiducial values a — 1, f3 — 1, = 1 
and Q = 0. As Ky increases, the growth rate o~ x becomes larger and the short wavelength 
perturbations become more and more unstable. The most unstable mode (the wavenumber 
at which a x attains its maximum value) increases as k\\ increases. We note that k XjCr i t , the 
critical wavenumber beyond which the system becomes stable is different for k\\ = and 
K\\ 7^ 0. This remind us about the effect of diffusion of CR on Jeans instability of uniform 
unperturbed state. When CR and gas are tightly bound together, the stability criterion 
depends on the sum of the gas and CR pressures. If there is CR diffusion (no matter how 
small it is), the stability criterion depends on the gas pressure only. The only effect of CR 
is on the growth rate and not the stability criterion. 

Fig. 2 (bottom) shows the dependence of the maximum growth rate (T Xiinax on k\\. o~ X:Uiax 
ranges from ~ 0.53 to ~ 0.65 in the range of < km < 200. When km < 10, o- x max increases 
rapidly with Ky, then it levels off to almost a constant when Ky > 100. 

Fig. 3 shows the dispersion relation of the case k x = for different km. The other 
parameters take the fiducial values a — 1, /3 — 1, — 1 and Q — 0. o~ y does not depend on 
Ky at all, which is understandable as CRs diffuse along magnetic field lines only. 

Fig. 4 shows the dispersion relation of general perturbations for Ky = 100 (left), 1 
(middle), and 0.1 (right). The other parameters take the fiducial values a — 1, (3 — 1, Zd — 1 
and Q = 0. Since o~ Xtmax decreases as km decreases while o- ymax is independent of km, it is 
expected that the most unstable mode shifts from k y = to k x = as Ky decreases. In fact, 
the transition occurs somewhere around Ky ~ 1 (Fig. 4 middle). As the coupling between 
CR and gas increases, the short wavelength perturbations in the unperturbed magnetic field 
direction are inhibited. 

Although the system is stable for large k y when k x = 0, it is unstable when k x is small 
but not zero. Moreover, the growth rate at large k y and small k x becomes smaller as Ky 
decreases. 



4.2. Dependence on Ratio of CR Pressure to Gas Pressure (3 

Fig. 5 shows the dispersion relation of the case k y = for different (3 (the ratio of the 
unperturbed CR pressure to gas pressure). The other parameters take the fiducial values 
a — 1, Ky = 100, Zd — 1 and fl = 0. a Xtmax increases (decreases) with (3 when j3 < 3 
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(> 3). The critical wavenumber fc XiC rit increases with (3 for (3 < 3, and it remains roughly the 
same for (3 > 3. As mentioned in previous subsection that CR should not affect the Jeans 
stability criterion of uniform unperturbed state when the diffusion coefficient is nonzero. 
The dependence on (3 of the critical wavenumber in this case comes from the dependence of 
the scale height of the unperturbed nonuniform state on the unperturbed CR pressure. 

Fig. 6 shows the dispersion relation of the case k x = for different f3. The other 
parameters take the fiducial values a — 1, k\\ = 100, z d — 1 and O = 0. a V;max decreases as 
(3 increases and levels off at (3 > 10. 

Fig. 7 shows the dispersion relation of general perturbations for (3 — 1 (left), 40 (middle), 
and 100 (right). The other parameters take the fiducial values oc — 1, k\\ = 100, z d — 1 and 
O = 0. At small (3 the most unstable mode has k y = 0. Around (3 ~ 40, a„ max and 

the most unstable mode shifts to k x = for larger (3. Moreover, the growth rate at large k y 
and small k x decreases slightly as (3 increases. 

4.3. Dependence on Half-Thickness z d 

Fig. 8 shows the dispersion relation of the case k y = for different z d (the half-thickness 
of the unperturbed gas disk). The other parameters take the fiducial values a — 1, f3 — 1, 
k\\ = 100 and O = 0. Both a XtIDax and k XtCT i t increase as z d increases. 

Fig. 9 shows the dispersion relation of the case k x = for different Zd- The other 
parameters take the fiducial values a — 1, km = 100, — 1 and O = 0. As Zd increases, 
°"?/,max gradually increases and /c^crit decreases quite a lot. 

Fig. 10 shows the dispersion relation of general perturbations for z d = 2 (left), 0.67 
(middle), and 0.35 (right). The other parameters take the fiducial values a = 1, j3 = 1, 
Kn = 100 and O = 0. When the z d is large (small) the most unstable mode occurs at ky — 
(k x = 0). When the perturbations develop into elongated fragments, they prefer to align 
perpendicular (parallel) to the unperturbed magnetic field. The transition between parallel 
and perpendicular fragmentation occurs around z d = 0.67. 

4.4. Dependence on Angular Velocity O 

Fig. 11 (left) shows the dependence of the maximum growth rate on (3 for different 
O (rotational angular velocity of the disk). The other parameters take the fiducial values 
a — 1, (3 — 1, k\\ — 100 and z d = 1. We note that the maximum growth rate is mostly 
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o"a;,max for the cases we studied here. As f3 increases, the maximum growth rate increases 
slightly until (3 ~ 3 for f2 = 0.0, (3 ~ 5 for VL = 0.5 and (3 ~ 9 for f2 = 1.0, then it decreases 
with (3. However, in the case of Q = 0.0, the maximum growth rate shifts to more or less 
constant at (3 > 40, because the maximum growth rate shifts from a x ^ max to <j yraaiX at (3 ~ 40. 
When we increase Q, the shift from a X)max to o~ y)inax occurs at larger and larger (3. When 
Q > 0.26 this shift of unstable modes never substantiates, and the k x = mode becomes 
stabilised (see Fig. 11, right). The growth of perturbations in this direction is stabilised by 
the Coriolis force. Thus when f2 is large enough the elongated filaments always align along 
the unperturbed magnetic field. 



5. SUMMARY AND DISCUSSION 

This is the first attempt to study how cosmic rays affect the Parker- Jeans instability of 
magnetized self-gravitating gaseous disks. We adopt a two-fluid model for the cosmic-ray- 
plasma system, in which the propagation of cosmic rays includes advection and diffusion. 
In this work we consider slab geometry with magnetic field lies along the slab, and cosmic 
rays diffuse only along the magnetic field lines. Moreover, we neglect external gravitational 
force. We perform standard linear perturbation analysis and work out the cases of symmetric 
boundary conditions at z = (Eq. (52), antisymmetric modes are less unstable). In general 
the most unstable mode occurs at either k y = or k x = 0, where the unperturbed initial 
magnetic field lies in the x-direction. Thus the perturbations will develop into elongated 
fragments with axes either perpendicular to (k y = 0) or parallel to (k x = 0) the initial 
magnetic field. 

The coupling between CRs and gas is contained in the diffusion coefficient in our model. 
In the case of k y = 0, the growth rate of the instability increases as k\\ increases (i.e., the 
coupling becomes weaker), and is levelled off at large «y (see Fig. 2). Recall that the thermal 
gas pressure is the stabilizing factor in both Jeans instability and Parker instability. One 
may interpret CR as a gas that can slip through the thermal gas according to the value of 
the diffusion coefficient or the strength of the coupling. During the process it exerts part 
of its pressure onto the thermal gas. In effect the thermal gas pressure is augmented and 
the Parker- Jeans instability becomes less unstable. We should point out that the instability 
cannot be stabilized by just decreasing the diffusion coefficient (it is less unstable only), 
except when the diffusion coefficient is exactly zero and the CRs and thermal gas are tightly 
bound together. See Fig. 2 (top) for detail. In the case of k x — 0, the dispersion relation 
does not depend on rey (see Fig. 3). This reflects the assumption that CRs diffuse along the 
magnetic field lines only. 
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Keeping other parameters fixed, the increase of j3 (the ratio of unperturbed CR pressure 
to gas pressure) has both stabilizing and destabilizing effect (see Fig. 5). The wavenum- 
ber beyond which the system is stable increases as (3 increases (it approaches a maximum 
value when /3 > 3), i.e., the unstable range increases. The maximum growth rate increases 
(slightly) with (3 when f3 < 3 but decreases with (3 when j3 > 3. We deem that the Parker 
mode is responsible for the increase in the unstable range and also for the increase in maxi- 
mum growth rate when (5 < 3, while the Jeans mode is responsible for the decrease in max- 
imum growth rate when (5 > 3. For comparison, Parker instability under constant external 
gravity (and no self-gravity is considered) is usually studied in an exponential atmosphere. 
As j3 increases, both unstable range and maximum growth rate increase until some maxi- 
mum values (see the appendix for a brief discussion). The effect of CRs on Parker instability 
under constant external gravity has, of course, been studied before, for instance, the original 
work by Parker (1966), and more recent works by Ryu et al. (2003) and Kuwabara et al. 
(2004) which have taken the diffusion of CR into account. However, one may notice that in 
Parker (1966) and Ryu et al. (2003), both unstable range and maximum growth rate of the 
instability increase without limit as (5 approaches oo. We should point out that the apparent 
increase without bound is just an artefact of normalisation. In these works the wavenumber 
and the growth rate are normalised to 1/H and u/H, where H is the scale height of the 
exponential atmosphere which is proportional to the total unperturbed pressure, and u is 
a characteristic speed (say, the isothermal sound speed). Thus the normalisation decreases 
as j3 increases. The actual unstable range and maximum growth rate are always finite (see 
appendix). 

Nagai et al. (1998) showed that the direction of the longitudinal axis of the elongated 
filaments formed by the growth of the perturbation with respect to the magnetic field depends 
on the half thickness of the disk or the external pressure. For thicker (thinner) disk the growth 
rate of perturbations in the x direction is larger (smaller) than those in the y direction. Here 
we demonstrate that the alignment can also be affected by CR diffusion. Fig. 12 shows the 
dependence of the growth rate on the half thickness of the disk Zd and the relation between 
o^max and <J y , max for different k\\ in the case of a — 1 and (5 — 1. As long as a x>max < cr^max 
the filaments will align along the magnetic field. As shown in the figure, when Zd < 0.65 
or Zd > 1.1, then 0" Xjmax is either always smaller or larger than <J y ^ nax and the diffusion of 
CR does not affect the direction of the filaments. However, when 0.65 < Zd < 1.1 (the grey 
region in the figure), the alignment depends on k\\. At a particular Zd the filaments prefer to 
align along the magnetic field when k\\ is small enough (i.e., the coupling is strong enough), 
and vice versa. 

The CR pressure also affect the alignment of filaments. Fig. 13 shows the dependence 
of the growth rates <7 x ,m&x and a y!iaax on (5 for three thickness Zd = 3, 1, and 0.66. Whenever 
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there is a crossover between the curves a Ximax and a VjUiax for the same Zd, the direction of the 
filaments depends on the CR pressure. In the case of Zd = 0.66, cr Xtmax is always smaller than 
<%max- For zj, = 1 and Zd — 3, <7 X)max < o^max when j3 > 40 and /3 > 550, respectively (also 
see Fig. 7 for the case of Zd — 1). In fact, for larger ^ crossover still exists for large enough 
(3, but the difference between <7 Xjmax and c^max is very small in those range of f3. Therefore, 
the alignment of the elongated filament depends on CR pressure when the thickness of the 
disk is large enough (z d > 0.66). 

The dependence of the eigenfunction 5p on several factors is shown in Fig. 14. The top 
panel shows the dependence on the CR diffusion coefficient k\\. The solid line shows Sp with 
K|| = 200 at the most unstable wave number (k x , k y ) = (0.5, 0), and the dotted line shows 
bp with any k\\ < 1 at the most unstable wave number (k x ,k y ) = (0, 0.35). The middle 
panel shows the dependence on the half thickness of the disk z d . The solid line shows the 
5p with Zd = 3 at the most unstable wave number (k x , k y ) = (0.59, 0) and the dotted line 
shows 5p with Zd — 0.1 at the most unstable wave number (k x , k y ) = (0, 2.94). The bottom 
panel shows the dependence on the ratio of CR pressure to gas pressure f3. The solid line 
shows the 5p with f3 — 3 at the most unstable wave number (k x ,k y ) = (0.51, 0) and the 
dotted line shows 5p with (3 = 1000 at the most unstable wave number {k x , k y ) = (0, 0.29). 
The solid line in each panel shows the case when the elongated filaments have the tendency 
to align their longitudinal axes perpendicular to the unperturbed magnetic field, while the 
dotted line in each panel shows the case when the elongated filaments prefer to align along 
the unperturbed magnetic field. Fig. 14 shows that the perturbation of density concentrate 
around the mid-plane if the most unstable mode is k x = (the dotted line), while it spread 
more evenly if the most unstable mode is k y = (the solid line). 

Numerical computations were carried out on the VPP5000 at the National Astronomical 
Observatory, Japan (NAOJ) under the projects wtk03b, and the SX-6 at the National insti- 
tute of Information and Communications Technology (NICT), Japan. CM. Ko is supported 
in part by the National Science Council of Taiwan, by grants NSC-92-2112-M-008-046 and 
NSC-93-2112-M-008-017. 



A. PARKER INSTABILITY IN AN EXPONENTIAL ATMOSPHERE 

In this appendix we briefly describe the dependence of the unstable range and the max- 
imum growth rate of the Parker instability under constant external gravity in an exponential 
atmosphere. 

In the unperturbed equilibrium state the external gravity g = — ge z and the magnetic 
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field B = £>e x . The hydrostatic equilibrium equation (g is constant) 

dP t 



dz 



-pg, 



(Al) 



gives an exponential atmosphere if P t oc p. Specifically, P t = P t0 exp[— z/H] and p = 
p exp[—z/H]. The scale height H = P t /(pg) = P t0 /(p g) is a constant. The total pressure 
P t includes thermal gas pressure, magnetic pressure and other pressures such as cosmic ray 
pressure. 

To study the linear stability of the system, one may perform standard normal mode 
analysis. We present the result of a particular simple case: only gas and magnetic field 
are actively involved in the dynamics of the system. The other components, e.g., cosmic 
rays, just provide the background pressure against the external gravity in the unperturbed 
equilibrium state. This is the case when cosmic ray diffusion coefficient is very large, i.e., the 
coupling is vanishing small. Let the perturbations be proportional to exp(<rt + ik x x) (forget 
about perturbations in y-direction to further simplify discussion) . The dispersion relation of 
the lowest mode in z-direction is 



- 4 + (V A 2 + C 2 ) (hi + j^) - 2 + klCl 



1 



AH 2 



V 2 + N 2 



0, 



(A2) 



where Va = B/y/p p and C s = a/ 7 g P g / P are the Alfven velocity and adiabatic sound speed, 
respectively; N = \Jg~j H — g 2 /C 2 is the Brunt- Vaisalaa frequency (with contributions from 
magnetic field and cosmic rays through H). The system is susceptible to Parker instability 
if 

AH 2 H 



N 2 + 



C 2 AH 2 



< 0. 



(A3) 



To facilitate discussion, we suppose P B = B 2 /2p = aP g , P c = f3P g and P g = pu 2 (the 
isothermal sound speed u is constant). Thus P t = P g +P B +P c = (l+a+(3)P g = (l+a+(3)pu 2 . 
A natural time scale of the system is u/g and a natural length scale is u 2 /g. In these 
units the parameters in Eq. (A2) become: V£ = 2a, C 2 = 7 g , H = (1 + a + fi) and 
A^ 2 = —(1 + a + j3 — 7 g )/(l + ol + /3)/7 g . The growth rate and wavenumber become oujg 
and k x u 2 /g. The instability criterion Eq. (A3) can then be written as (cf., Parker 1966) 



;i + a + (3) 2 > 7 g (l + \a + 0) 



(A4) 



In these units it can readily be shown that the critical wavenumber beyond which the system 
is stable is given by 



u 



1 



2«7g 



lg {l + \a + l3) 
(1 + a + p) 2 



(A5) 
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and the maximum growth rate is given by the smaller root of 



(2a - 7 g ) (J^max- -2(2a + 7 g ) 1 




7g (l + iq + /?) / 
(l + a + (3) 2 V 




+ 



7g (l + fq + /3) 
(l + a + /3) 2 



0. 



(A6) 



In fact, {k x ^ cr i t u 2 / g) 2 and (o" x ,max M /fi') 2 increase monotonically with /3. Moreover, as /3 
(k x , crit u 2 / g) 2 -> l/2«7 g and (cr x ^u/ g) 2 -> l/(v / 2a + ^/Tg). 
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Fig. 1. — Left: Gas pressure in the unperturbed equilibrium state as a function of z (vertical 
coordinate) in the case of a = 1 and (5 — 1. Right Gravity in the unperturbed equilibrium 
state as a function of z. 

Fig. 2. — Top: The dependence of the dispersion relation at k y = on k», the CR diffusion 
coefficient. a x is the growth rate when k y = and k x is the wavenumber in the direction of 
the unperturbed magnetic field. The other parameters are a — 1, f3 — 1, Zd — 1 and Q = 0. 
Bottom: The dependence of the maximum growth rate cr^max on k\\. 

Fig. 3. — The dependence of the dispersion relation at k x = on k\\. a y is the growth rate 
when k x — 0. The other parameters are the same as in Fig. 2. 

Fig. 4. — The general dispersion relation for the Parker- Jeans instability against the effect 
of CR diffusion. Left: k\\ = 100, Middle: k\\ = 1, and iZig/ifc k\\ = 0.1. The other parameters 
are the same as in Fig. 2. 

Fig. 5. — The dependence of the dispersion relation at k y = on /3, the ratio of the un- 
perturbed CR pressure to gas pressure. o~ x is the growth rate when k y = 0. The other 
parameters are a — 1, k\\ — 100, Zd — 1 and Q = 0. 

Fig. 6. — The dependence of the dispersion relation at k x = on (5. a y is the growth rate 
when k x — 0. The other parameters are the same as in Fig. 5. 

Fig. 7. — The general dispersion relation for the Parker- Jeans instability against the effect 
of CR pressure. Left: (3 = 1, Middle: (3 = 40, and Right: (3 = 100. The other parameters 
are the same as in Fig. 5. 

Fig. 8. — The dependence of the dispersion relation at k y = on Zd, the half-thickness of the 
unperturbed gas disk (can be viewed as a measure of external pressure). o~ x is the growth 
rate when k y = 0. The other parameters are a — 1, (3 — 1, K\\ — 100 and Q = 0. 

Fig. 9. — The dependence of the dispersion relation at k x — on z d . a y is the growth rate 
when k x — 0. The other parameters are the same as in Fig. 8. 

Fig. 10. — The general dispersion relation for the Parker- Jeans instability against the effect 
of the half-thickness of the gas disk. Left: z d = 2, Middle: z d = 0.67, and Right: z d = 0.35. 
The other parameters are the same as in Fig. 8. 
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Fig. 11. — Left: The dependence of the maximum growth rate <7 max on (3 (the ratio of 
unperturbed CR pressure to gas pressure) at different angular velocity of the gas disk Q. 
The other parameters are a — 1, f3 — 1, k» — 100 and z d = 1. Right: The general dispersion 
relation in the case of a — 1, f3 — 1, kh — 100, z d — \ and Q = 0.5 



Fig. 12. — The dependence of the maximum growth rate on the half-thickness z d . cr x iaax is 
the maximum growth rate at k y = at a given k\\. <7 y , max is the maximum growth rate at 
k x = (a y<max is independent of k\\, see Fig. 3). The other parameters are a — 1, (5 — 1 
and = 0. In the grey region, CR diffusion determines which mode is more unstable, thus 
determines the alignment of the perturbations with respect to the unperturbed magnetic 
field. 



Fig. 13. — The dependence of the maximum growth rate on the ratio of CR pressure to gas 
pressure f3 at different z d . The other parameters are a — 1, k\\ = 100 and f2 = 0. 



Fig. 14. — The eigenfunction dp of the most unstable mode. Top: The dependence on the CR 
diffusion coefficient ku. The solid line is the case of kh = 200 at the wavenumber (k x , k y ) = 
(0.5, 0), and the dotted line is the case of any k\\ < 1 at the wavenumber (k x , k y ) = (0, 0.35). 
The other parameters are a = 1, (3 = 1, z d = 1 and f2 = 0. Middle: The dependence on 
the half-thickness z d . The solid line is the case of z d = 3 at (k x ,k y ) = (0.59, 0), and the 
dotted line is the of z d = 0.1 at (k x , k y ) = (0, 2.94). The other parameters are a = 1, f3 — 1, 
Ky = 100 and Q, — 0. Bottom: The dependence on the ratio of CR pressure to gas pressure 
(3. The solid line is the case of /3 — 3 at (k x , k y ) = (0.51, 0), and the dotted line shows the 
case of P — 1000 at (k x , k y ) = (0, 0.29). The other parameters are a — 1, k\\ = 100, z d = 1 
and Q = 0. 
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